#include "mt3d_cuda.h"

void mt3d::read_model(std::string filename)
{
    std::clog << "Reading Tetgen model ... ";
    gctl::read_Tetgen_node(filename, mesh_nodes_, gctl::Packed, &node_markers_);
    gctl::read_Tetgen_element(filename, mesh_elems_, mesh_nodes_, gctl::Packed, &elem_markers_);
    gctl::read_Tetgen_face(filename, mesh_faces_, mesh_nodes_, gctl::Packed, &face_markers_);

    elem_num_ = mesh_elems_.size();
    node_num_ = mesh_nodes_.size();
    face_num_ = mesh_faces_.size();

    // Sort out edge lists
    tet_edge tmp_edge;
    std::string tmp_str;
    for (size_t i = 0; i < elem_num_; i++)
    {
        for (size_t j = 0; j < 6; j++)
        {
            if (mesh_elems_[i].vert[le_index[2*j]]->id < mesh_elems_[i].vert[le_index[2*j+1]]->id) // edge index is from small to big
            {
                tmp_edge.set(mesh_elems_[i].vert[le_index[2*j]], mesh_elems_[i].vert[le_index[2*j+1]]);
            }
            else
            {
                tmp_edge.set(mesh_elems_[i].vert[le_index[2*j+1]], mesh_elems_[i].vert[le_index[2*j]]);
            }

            // Try to find the edge to see if it is already existed
            tmp_str = std::to_string(tmp_edge.vert[0]->id) + " " + std::to_string(tmp_edge.vert[1]->id);
            edge_index_iter_ = edge_index_map_.find(tmp_str);

            if (edge_index_iter_ != edge_index_map_.end()) // already existed
            {
                mesh_edges_[edge_index_iter_->second].host.push_back(mesh_elems_.get(i));
            }
            else // do not existed
            {
                tmp_edge.id = mesh_edges_.size(); // assgin edge index
                edge_index_map_[tmp_str] = tmp_edge.id; // add to the map
                mesh_edges_.push_back(tmp_edge); // add the edge to list
                mesh_edges_.back().host.push_back(mesh_elems_.get(i)); // add the first tetrahedron
            }
        }
    }

    edge_num_ = mesh_edges_.size();

    // Sort element to global edge index
    bool reversed;
    gctl::tetrahedron *tet_ptr;
    elem_edge_index_.resize(elem_num_, 6);
    elem_edge_reversed_.resize(elem_num_, 6);
    for (size_t i = 0; i < edge_num_; i++)
    {
        for (size_t j = 0; j < mesh_edges_[i].host.size(); j++)
        {
            tet_ptr = mesh_edges_[i].host[j];
            for (size_t k = 0; k < 6; k++)
            {
                if (tet_ptr->vert[le_index[2*k]]->id < tet_ptr->vert[le_index[2*k+1]]->id) // edge index is from small to big
                {
                    tmp_edge.set(tet_ptr->vert[le_index[2*k]], tet_ptr->vert[le_index[2*k+1]]);
                    reversed = false;
                }
                else
                {
                    tmp_edge.set(tet_ptr->vert[le_index[2*k+1]], tet_ptr->vert[le_index[2*k]]);
                    reversed = true;
                }

                if (tmp_edge == mesh_edges_[i])
                {
                    elem_edge_index_[tet_ptr->id][k] = i;
                    elem_edge_reversed_[tet_ptr->id][k] = reversed;
                    break;
                }
            }
        }
    }

    // Sort out boundary edges
    edge_on_boundary_.resize(edge_num_, false);
    for (size_t i = 0; i < mesh_faces_.size(); i++)
    {
        if (face_markers_[i] == bound_marker_) // boundary faces contain boundary edges
        {
            for (size_t j = 0; j < 3; j++)
            {
                if (mesh_faces_[i].vert[j]->id < mesh_faces_[i].vert[(j+1)%3]->id) // edge index is from small to big
                {
                    tmp_edge.set(mesh_faces_[i].vert[j], mesh_faces_[i].vert[(j+1)%3]);
                }
                else
                {
                    tmp_edge.set(mesh_faces_[i].vert[(j+1)%3], mesh_faces_[i].vert[j]);
                }

                // Find the edge
                tmp_str = std::to_string(tmp_edge.vert[0]->id) + " " + std::to_string(tmp_edge.vert[1]->id);
                edge_index_iter_ = edge_index_map_.find(tmp_str);

                if (edge_index_iter_ != edge_index_map_.end()) // must to be existed
                {
                    edge_on_boundary_[edge_index_iter_->second] = true;
                }
                else throw gctl::runtime_error("Invalid edge found in .face file.");
            }
        }
    }

    std::clog << "done\n";
    std::clog << node_num_ << " nodes, " << elem_num_ << " elements, " << edge_num_ << " edges.\n";
    return;
}
